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ABSTRACT 

The dynamic spectrum of a radio pulsar is an in-line digital hologram of the ionised interstellar 
medium. It has previously been demonstrated that such holograms permit image reconstruc- 
tion, in the sense that one can determine an approximation to the complex electric field values 
as a function of Doppler-shift and delay, but to date the quahty of the reconstructions has 
been poor Here we report a substantial improvement in the method which we have achieved 
by simultaneous optimisation of the thousands of coefficients that describe the electric field. 
For our test spectrum of PSR B0834+06 we find that the model provides an accurate represen- 
tation of the data over the full 63 dB dynamic range of the observations: residual differences 
between model and data are noise-like. The advent of interstellar holography enables detailed 
quantitative investigation of the interstellar radio-wave propagation paths for a given pulsar at 
each epoch of observation; we illustrate this using our test data which show the scattering ma- 
terial to be structured and highly anisotropic. The temporal response of the medium exhibits 
a scattering tail out to beyond 100 /is and a pulse arrival time measurement at this frequency 
and this epoch of observation would be affected by a mean delay of 15 /xs due to multipath 
propagation in the interstellar medium. 

Key words: techniques: interferometric — pulsars: general — pulsars: individual: B0834-H06 
— scattering — ISM: structure — turbulence 



1 INTRODUCTION 

With modern instrumentation pulsar dynamic spectra can be 
recorded at high spectral and temporal resolution yielding a dataset 
with a large information content from just one observation. For ex- 
ample: if we observe a pulsar for one hour, sampling the spectrum 
with 1 kHz channels every 10 seconds, over a total bandwidth of 
100 MHz, then we have ~ 4 x 10^ independent flux measurements. 
And if the pulsar is bright and the telescope is large then each of 
these measurements can have signal-to-noise ratio of order unity 
implying a total information content of potentially ~40 Mbit. This 
information relates primarily to multipath scattering of the radio- 
waves by the ionised Interstellar Medium (ISM), as it is this scat- 
tering which gives rise to the observed interference fringes. It may 
be that these paths are determined by random irregularities - e.g. 
caused by turbulence - in the ISM. In that case any given dynamic 
spectrum contains information about those random irregularities, 
and there is no strong motivation to study the individual spectra in 
detail as they reflect particular realisations of a stochastic process. 
But some pulsar dynamic spectra exhibit a high level of order in 
their fringe patterns — see Rickett (I99I) for an overview. This 



fact has been emphasised by consideration of the two-dimensional 
power spectra of the dynamic spectra, wherein power is often seen 
to be concentrated in parabolic arcs and inverted arclets (Stinebring 
et al 2001). There is a consensus that this phenomenon arises di- 
rectly from the geometry of the scattering process (Stinebring et al 
2001; Cordes et al 2006; Walker et al 2004), with waves scattered 
through an angle 9 experiencing a Doppler-shift proportional to one 
component of and a delay proportional to • 9. In cases where the 
parabolae are very sharp it has been argued that the scattering is 
highly anisotropic (Walker et al 2004; Cordes et al 2006). Sharply 
defined arcs/arclets also require that the scattering arises in a region 
of small extent along the line-of-sight (Stinebring et al 2001), so it 
is not distributed turbulence but discrete, localised structures which 
are responsible for these features. However, the physical nature of 
the scattering media remains obscure. Consequently there is now 
some incentive to explore the information content in individual dy- 
namic spectra, in order to build a detailed picture of the scattering 
structures. Further motivation for investigating individual dynamic 
spectra is provided by studies of the pulsars themselves: if the prop- 
erties of the scattering screen are accurately known it is possible to 
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use the screen for very high resolution imaging of tiie pulsar mag- 
netosphere (Wolszczan and Cordes 1987; Gwitin et al 1997; Walker 
and Stinebring 2005, WS05 henceforth). Precision pulsar timing 
programs (e.g. Manchester 2007) also provide an incentive for un- 
derstanding the particular interstellar propagation paths which con- 
tribute to individual observations — if the propagation delays re- 
main uncorrected in the data they constitute a potentially large sys- 
tematic error in pulse arrival time measurements. 

It has previously been demonstrated that one can itcratively 
construct a model of the electric field as a function of radio- 
frequency and time, starting from the observed dynamic spectrum 
(WS05). The procedure for achieving this is largely equivalent to 
determining a phase for the electric field in each pixel of the dy- 
nanuc spectrum, because a noisy estimate of the field amplitude 
is given directly by forming the square root of the observed in- 
tensity. This situation is common to the broad class of problems 
known as "phase-retrieval problems", which are well known in the 
optics literature (e.g. Fienup 1982). However, the method of solu- 
tion demonsttated for pulsar dynamic spectta appears to be new; 
it exploits sparseness of the power distribution in the Fourier do- 
main, and a solution is obtained iteratively by adding discrete new 
field components in such a way as to minimise the differences be- 
tween the model dynamic spectrum and the data. Conceptually the 
process has a strong similarity to the CLEAN algorithm (Hogbom 
1974) which is commonly used in radio astronomical imaging for 
deconvolving the synthesised telescope beam from a "dirty" image 
of the sky; CLEAN usually works well if the image is only sparsely 
populated with errussion. The connection between the two algo- 
rithms is reinforced when we recognise that determining the elec- 
tric field from the dynamic spectrum is also equivalent to a decon- 
volution. The dynamic spectrum as a function of radio-frequency 
and time is simply 7(^, t) = U* (v, i) U {v, t), where U is the elec- 
tric field;_in the^Fourier domain this relationship becomes a con- 
volution I = U* ^ U and so we are deconvolving the Fourier 
ttansform of the electric field from its complex conjugate. 

Some comments about nomenclature are appropriate at this 
point. Any process which allows us to reconstruct the electric field 
which gave rise to a recorded Mnge pattern is sensibly termed 
"holography", and the recorded fringe pattern (i.e. the dynamic 
spectrum in our case) is the hologram. In the case considered here 
it is "digital holography" because the reconstruction is done in soft- 
ware, and "in-line" because the object (i.e. the scattering screen in 
our case) is ttansparent and sits in the beam which forms the refer- 
ence wave. In-line holography is sometimes called "Gabor hologra- 
phy" because it is the arrangement which was originally conceived 
by the inventor of holography, Dennis Gabor. 

To date the fideUty of electric field reconstructions from pulsar 
dynamic spectra has been poor. Although the model of WS05 suc- 
cessfully reproduces the general appearance of their test dynamic 
spectrum there are substantial quantitative differences between the 
model and the data. These errors are seen when the model "sec- 
ondary spectrum" (i.e the power spectrum of the dynamic spec- 
trum) is compared to the data: the data exhibit a dynamic range of 
63 dB, but the difference between the WS05 model secondary spec- 
trum and the secondary spectrum of the data has a dynarruc range 
of 47 dB, implying that the reconstruction has correctly captured 
only the top 16 dB of the data. Here we report modifications to the 
reconstruction process which have yielded dramatic improvements 
to the accuracy of the electric field model; on the same test dataset 
as used by WS05 we find that our improved model captures the full 
dynamic range of the data. These gains were achieved by simulta- 
neous optimisation of the thousands of parameters describing the 



wave interference phenomenon, and by simultaneous optimisation 
of the hundreds of parameters which describe the intrinsic temporal 
flux variations of the pulsar. 

This paper is organised as follows: in the next section we detail 
the improvements we have made to the reconstruction algorithm 
described by WS05; in §3 we present the results we have obtained, 
using the same test data employed by WS05; and in §4 we con- 
sider what our test data tell us about the interstellar medium, with 
emphasis on the temporal response of the scattering medium. 



2 OPTIMISATION OF THE E-FIELD MODEL 

In WS05 we described an algorithm for modelling the electric field 
structure in the Fourier domain conjugate to the dynamic spec- 
trum. The latter is recorded as a function of radio-frequency, v, 
and time, t, and the corresponding conjugate variables are delay, t, 
and Dopplcr-shift, ^. The WS05 algorithm proceeds from a grid of 
noisy measurements of the electric field envelope, \U{v,t)\'^, to a 
model of U{t, oj) in terms of discrete wave components, j: 

U{T,U!)=^UjS(T-Tj)S(w-Wj), (1) 
j 

and the components (characterised by Tj , cjj and the complex num- 
ber Uj) are chosen one-by-one so as to minimise the differences 
between model and data. With ~ 9000 components the model re- 
ported by WS05 gives a good visual match to the data, but the 
residuals are large in comparison with the noise level, implying 
large systematic errors. An inspection of the differences between 
the observed secondary spectrum (i.e. the power spectrum of the 
dynamic spectrum) and the model secondary spectrum - the two 
quantities shown in the lower panel of figure 1 in WS05 - reveals 
that the discrepancies occur predominantly at the same locations in 
delay-Doppler space where there is already power present in both 
model and data. This suggests that the systematic errors introduced 
by the WS05 algorithm are not due to incorrectiy placed compo- 
nents (i.e. wrong Tj,u}j), or an insufficient number of components 
in the model, but rather due to errors in determining the various Uj . 

It was noted in WS05 that global optimisation of the {uj} is 
desirable, in order to reduce systematic errors in the model, but dif- 
ficult to achieve because of the large number of free parameters 
which would be involved in such an optimisation. In particular in- 
version of a 10* X 10* non-sparse matrix - which is perhaps the 
most obvious approach to solution of the linearised least-squares 
optimisation - is computationally challenging. Furthermore any ap- 
proach based on solution of the linearised problem would require 
iteration in order to solve the full non-linear optimisation problem. 
Fortunately there are easier methods - see, for example, Nocedal 
and Wright (1999) - and we have used one of these to optimise the 
WS05 electtic field model as we now describe. 

The method which we employed is a quasi-Newton method, 
in which a demerit function S is minimised by seeking succes- 
sively closer approximations to the solution of dS/dxi = for all 
parameters Xi over which we wish to optimise. Newton's method 
requires knowledge of the Hessian (i.e. all the second derivatives 
d'^ S / dxidxj), which is computationally expensive when a large 
number, iV, of parameters is involved. By contrast, quasi-Newton 
methods proceed by approximating the Hessian; information from 
current and previous iterations is used to update the approxi- 
mate Hessian for subsequent iterations yielding, in effect, a finite- 
difference representation of the local curvature of S. The update 
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scheme which we used is the BFGS (Broyden-Fletcher-Goldfarb- 
Shanno) update. Specifically: we used the L-BFGS-B algorithm 
(Byrd et al 1995), which is a "Limited memory" implementation of 
BFGS, using line-search minimisation, in which Bounds are per- 
mitted on the parameters Xi. The term "Limited memory" here 
refers to the fact that the N x N Hessian is replaced by the outer 
product of two m X N matrices, with m being the number of prior 
iterations which are employed in constructing the update. Because 
m ~ a few, and m does not grow with A*', the storage requirements 
of the algorithm are only modest and grow linearly with A'^. The 
L-BFGS-B code was written by specialists in the field of numerical 
optimisation; it is freely available as a set of FORTRAN subrou- 
tine.']^ In order to make use of this code it is necessary for the user 
to supply routines which evaluate the demerit function, S, and its 
partial derivatives, dS/dxi, with respect to all the parameters, Xi, 
over which we wish to optimise. Given these inputs the L-BFGS-B 
code will search for a minimum in S. If a minimum is found by 
L-BFGS-B it is not guaranteed to be a global minimum. The L- 
BFGS-B package is well documented; there are clear instructions 
on how to use the code and simple drivers are included in the pack- 
age so that the user can verify correct performance on their own 
computer. For our application we used a model with two sets of 
parameters: the various Uj, each of which is described by two un- 
bounded real variables, representing the real and imaginary parts, 
and a set of positive-definite real numbers describing the intrinsic 
pulsar flux as a function of time, {fk}- The inclusion of the various 
parameters fk : J? V fc demands that the optimisation soft- 
ware be able to handle variables which are bounded, as is the case 
for the L-BFGS-B package. The WS05 algorithm does not attempt 
to solve for the fk explicitly but simply applies a Fourier-domain 
filter to the data in an attempt to remove the intrinsic pulsar flux 
variations. Because there is no clear-cut distinction between intrin- 
sic flux variations and those due to wave interference, the procedure 
used by WS05 is quite crude and in the present work we use it only 
as a starting point (as per item (i) in §3.1 of WS05). 

Our first attempt at improving on the WS05 approach was to 
take the output of the WS05 algorithm and use it as the starting 
point of an optimisation with L-BFGS-B. The results were good, 
yielding a model which captured much more of the dynamic range 
in the data, and had a lower value of the reduced statistic. This 
result was encouraging, but it was clear that we could do better: the 
components which are identified at each iteration in the WS05 al- 
gorithm depend on the electric field model at that point, so the sys- 
tematic errors in the WS05 model are not completely eliminated by 
post facto optimisation — spurious components remain in the opti- 
mised model, albeit at a low level. If, on the other hand, the electric 
field model at each iteration of the WS05 algorithm is optimised 
(using L-BFGS-B) then these spurious features can be minimised. 
This approach has the additional benefit of simplicity in the pro- 
cessing of the data, requiring only one pass through the data and 
one set of code. We therefore implemented a new iterative decom- 
position algorithm, in FORTRAN, which employs the L-BFGS-B 
package to optimise the model Uj and fk . The new algorithm dif- 
fers from WS05 in the following ways: 

(i) The parameters fk (describing the intrinsic pulsar fluxes) 
are optimised once for every 100 new field components which are 
picked. The optimisation over {fk} is done separately from that 
over {uj}. 
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(ii) The number of new electric field components picked at each 
iteration is given by the integer part of 1 + 7Vc/100, where Nc is the 
current number of electric field components. So initially only one 
new component is picked at each iteration, and when the model 
contains a large number of components the fractional increase per 
iteration is 1%. 

(iii) When the number of electric field components exceeds 100, 
the algorithm is free to pick components which have t < 0. 

(iv) The algorithm terminates when the reduced (i-e. P^i' 
degree of freedom) reaches unity, or when the reduced reaches 
a minimum — whichever occurs first. 

Because of the high dynamic range of the data it is important to 
maintain high precision in the optimisation, so we set the L-BFGS- 
B parameter named "factr" (which measures the precision in units 
of machine precision) to the value 10. The number of previous steps 
used in forming each BFGS update is m = 10. 



3 RESULTS 

A direct comparison of the results of WS05 with the new algorithm 
is possible by using the same test data as WS05 employed. Those 
data are shown in figure 1, along with the model generated by the 
new algorithm; both model and data are shown in the form of a 
dynamic spectrum and its power spectrum (the "secondary spec- 
trum"). The model was generated using the algorithm described in 
§2. Figure 1 also shows the residuals between data and model dy- 
namic spectra, and the power spectrum of those residuals; the resid- 
uals appear noise-like in both panels. Some quantitative measures 
of the success of the new algorithm are appropriate. The new algo- 
rithm achieves a reduced chi-squared value of Xr = 1-00 (this was 
the stopping criterion which was reached first), versus Xr ~ 1-19 
achieved by WS05, and it does so with only 8,000 electric field 
components versus 8,720 in WS05. Note, however, that the new 
algorithm does employ an additional 270 real numbers - one for 
each time sample - to describe the intrinsic pulsar flux variation 
with time; these numbers were in effect fixed in the WS05 algo- 
rithm by a Fourier-plane filter acting on the input data. Subtract- 
ing the model dynamic spectrum from the data, and then forming 
the power spectrum of the residuals gives a sensitive test of the 
fidelity of the model because it picks out correlated errors in the 
dynamic spectrum model. For the new algorithm the largest value 
of the residual power is only 1 1 dB above the mean noise power in 
the data, whereas the peak power in both data and model is 63 dB 
above the mean noise power in the data. Thus the algorithm is cer- 
tainly free of systematic errors over a range of 52 dB. 

In fact the new algorithm has achieved noise-limited perfor- 
mance and is capturing the full dynamic range of the data; to see 
this we need only examine the statistics of the noise power, shown 
as a histogram in figure 2. If the residuals were purely noise-like 
then the expected probability distribution function would be an ex- 
ponential, because the residual power is the sum of the squares of 
two independent variables each of which has the same Gaussian 
distribution. We can see from figure 2 that the residuals conform 
closely to this expectation; we can also see that the peak residual is 
not introduced by any systematic error in the modelling but rather 
it is simply the tail of the noise distribution. The residuals actually 
exhibit a slightly lower noise level than the data (dashed line); this 
can be understood by reference to figure 1. The noise in the data 
is estimated from the floor power level in the secondary spectrum, 
specifically an area in the corner of the secondary spectrum (away 
from any obvious signal power) is chosen and the mean power over 
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Figure 1. Data (top), model (middle), and residuals between data and model (bottom), for an observation of PSR B0834+06 in a 1.56 MHz band centred on 
321.00 MHz. The data were taken with the Arecibo radio telescope in conjunction with the WAPP backend signal processing units on MJD 53009; there are 
1024 spectral channels, and 270 time samples each of 10 seconds duration. The left-hand panels show dynamic spectra, while the right-hand panels show the 
corresponding secondary spectra (power spectra of the dynamic spectra); the range of the secondary spectra is ±50 mHz in Doppler-shift and ±327.6 /us in 
delay. Inverse grey-scale (black is peak intensity) is used in all cases; the transfer function is linear for the dynamic spectra, and logarithinic for the secondary 
spectra. All secondary spectra shown here have the same transfer function; the transfer functions for dynamic spectra are identical in the case of data and 
model, whereas the output range is stretched by a factor of five for the dynamic spectrum residuals in order to reveal their structure. In comparison with figure 
1 of WS05 note that the present figure includes the modulating effects of the intrinsic pulsar flux variations. We have chosen to display the full secondary 
spectra, including negative delays, even though these spectra are symmetric through the origin, so that the structure near zero delay can be better seen. 
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Power / Peak Power 

Figure 2. Histogram of the number of pixels per bin (of width 6 x 10~®) as 
a function of power, for the two-dimensional power spectrum of the resid- 
ual between the model dynamic spectrum and the data (solid line). Also 
shown is the corresponding histogram for the data (dashed line). The resid- 
ual power distribution corresponds very closely to the expected exponential 
noise distribution. The mean power in the residuals is slightly less than that 
calculated from the data because the latter estimate includes a small contri- 
bution from signals which are incorporated into the model. 



this area is computed. Tlie model secondary spectrum also extiibits 
a floor power level, even though the model is intended to represent 
"signal" rather than "noise"; this power is not present in the resid- 
ual dynamic spectrum so the mean power level in the residuals is 
slightly lower than in the data. 

The model electric field strength determined by our new algo- 
rithm is shown in figure 3. Visually this result is similar to that ob- 
tained by WS05 (their figure 2), with the notable exception that in 
WS05 all field amplitudes were fixed at zero in the region r < (so 
that region is not displayed in their figure 2). On physical grounds 
this area is expected to be free of astronomical signals. However, 
in-Une holography generally involves a certain amount of confu- 
sion between the image and its conjugate, because the observable 
quantity is usually the intensity U*U which is invariant under the 
replacement U — > ?7*. In the present context the conjugate image 
appears at negative delays, because forming the complex conjugate 
of cxp[2Tri(i^Tj + iOjt] is equivalent to making the replacements 
Tj —Tj, t^j —> —^j- In our figure 3 it can be seen that there is 
little trace of an inverted parabola in the lower half of the figure - 
only weak components can be seen in the region t < - indicat- 
ing minor confusion between the image and its conjugate (see also 
§4.1). This clean separation is another indication of the low level of 
systematic errors inherent in the new algorithm. In the same vein 
we note that there is little power evident in figure 3 near the r = 
locus (a horizontal Une tangent to the apex of the parabola), indicat- 
ing that the intrinsic pulsar flux variations, described by the various 
fk, have been accurately modelled. 

Although the origin of coordinates (i.e. r = 0, w = 0) in fig- 
ure 3 is in principle unknown (WS05), in the case of this particular 
dataset it appears sensible to assign the origin to the centre of the 
image because this is where the largest amplitude field component 
is located and the great majority of the intensity in figure 3 lies 
above this point — consistent with the brightest image component 
being an unscattered wave. 

Having now reached the point where we can form models 
which are a good match to the data it is appropriate to draw some in- 




Figure 3. Model electric field amplitudes corresponding to the model spec- 
tra shown in figure 1. Here the amplitudes are shown in grey-scale, as a 
function of Doppler-shift and delay, with a logarithmic transfer function; 
peak intensity is black. The Doppler-shift ranges over a total of 100 mHz 
(270 pixels) and the delay ranges over a total of 655 (1024 pixels). Ab- 
solute values of delay and Doppler-shift are in principle unknown but it ap- 
pears sensible to assign the origin of coordinates to the centre of the image 
(which is the peak field amplitude) in this particular case. It is evident from 
this figure that the image has been well separated from its complex conju- 
gate, which would appear as an inverted parabola. Also notable is the low 
level of power around zero delay, which indicates that the intrinsic pulsar 
flux modulations have been accurately determined by the algorithm. 

ferences about the properties of the scattering medium which gives 
rise to the data shown in figure 1. Excepting the temporal analysis 
in §4.1 our discussion is only quaUtative; that is because we have 
constructed our holographic image in delay-Doppler space (i.e. the 
Fourier space conjugate to the frequency-time space in which the 
data are recorded), whereas progress in understanding the scatter- 
ing medium relies on a knowledge of the image in two-dimensional 
spatial (angular) coordinates and there is no simple way of proceed- 
ing from the former to the latter. 



4 PROPERTIES OF THE SCATTERING MEDIUM 

An important qualitative aspect of figure 3 is that power in the 
model electric field is tightly concentrated around a parabolic lo- 
cus, with delay proportional to the square of the Doppler-shift; this 
much was already evident in WS05. This confirms that the under- 
lying scattering is highly anisotropic, with a scattered image which 
is very much longer than it is wide — a conclusion which has pre- 
viously been arrived at from consideration of the properties of ob- 
served secondary spectra for this and other pulsars (Walker et al 
2004; Cordes et al 2006; Trang and Rickett 2007). 

We can also see from figure 3 that some parts of the parabola 
show high intensity levels while others are almost completely de- 
void of signal, and in some places the high-low transitions are fairly 
abrupt. Abrupt changes are suggestive of weU-deflned boundaries 
to the scattering regions. Four intensity concentrations are seen at 
large delay and positive Doppler-shift; these correspond to the fea- 
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(cf. equations 1 and 2 of WS05). The corresponding field intensity 



100 

Delay (/^s) 

Figure 4. The mean intensity impulse response of the scattering medium: 
top panel, as determined from the holographic image shown in figure 3; 
bottom panel, as determined from the holographic image shown in figure 3 
with low amplitude coefficients (\uk \ < 0.004) set to zero. 



tures named "A,B,C,D" by Hill et al (2005) in a secondary spectrum 
analysis of data spanning more than three weeks. These features 
are presumably due to localised plasma concentrations; it is not yet 
clear whether these concentrations are diffracting or refracting the 
radio-waves into the telescope. Feature "A" lies significantly above 
the locus of the parabola; this extra delay could be the wave-speed 
(dispersive) delay of a high column-density structure — an inter- 
pretation which can in future be tested by comparing data at two 
different frequencies obtained at the same epoch. 



4.1 Temporal response of the medium 

Inferring the spatial structure of the scattering medium from our 
image of U is not a simple exercise and it is beyond the scope of 
this paper to attempt an in-depth analysis. On the other hand the 
holographic image shown in figure 3 is well suited to determin- 
ing the temporal response of the medium which is introduced by 
multi-path propagation. The electric field amplitude as a function 
of delay, t, and observing time, t, can be determined by forming 
the inverse Fourier transform of equation 1 with respect to ui: 



I{T,t)=U*U =A{t) + B{T,t), 



(3) 



where 

B{T,t) 



j) 5{t — Tk) exp [2m{ujk — Wj)*] : 

(4) 

describes the beating between waves of differing Doppler-shifts 
(but the same delay) and 



do; f/(r, w) 



(5) 



is the mean intensity impulse response function of the scattering 
medium. The function Aij) convolved with the intrinsic pulse pro- 
file yields (up to a normalisation factor) the average pulse profile 
for this observation. For our image, U, the mean intensity impulse 
response, A, is shown in the top panel of figure 4. For those delays 
where U includes one wave component which has a much larger 
amplitude than the other components at that delay the beat terms 
will all be relatively small and \B\ will be small compared with A. 
In general, however, B is not negligible; we defer consideration of 
B to later in this section. 

There are several recognisable features of the A{t) deter- 
mined from our model: the peak at zero delay is in accord with the 
physical expectation that a bright image should form at the delay 
minimum; the parabolic arc seen in figure 3 is responsible for the 
extended scattering tail stretching out beyond r = 100 /is; and the 
peaks near t = 140, 160, 230, 270 /is in figure 4 correspond to the 
features labelled "A, B, C, D", respectively, by Hill et al (2005) — 
these features are seen as discrete intensity concentrations at large 
delay and positive Doppler-shift in figure 3. 

At large negative delays in the top panel of figure 4 wc sec an 
intensity level ~ 10~^, whereas physically we expect zero inten- 
sity because wave propagation can only introduce positive group- 
delays. This is simply due to noise in the reconstruction; note also 
that a similar floor intensity level is present at large positive delays. 
The model electric field inevitably includes noise because there is 
no clear distinction between noise and signal contributions to the 
input data. By setting the low amplitude coefficients of U to zero 
we are able to eliminate much of this intensity floor. Specifically, 
setting the model coefficients Uj to zero for \ uj\ < 0.004 preserves 
all of the recognisable signal features in A(t) but largely removes 
the noise floor; the result is shown in the lower panel of figure 4. 

At small negative delays the intensity level in figure 4 is up 
to an order of magnitude higher than the noise floor. We interpret 
this as evidence of low-level confusion between the holographic 
image, U , and its complex conjugate, U* — as noted in §3 we 
expect this confusion to be present at some level. For applications 
where suppression of the conjugate image is critical it is possible 
to undertake the holographic image reconstruction entirely in the 
positive-delay half-space, but we note that this approach is expected 
to be problematic if the brightest image is not the image with the 
least delay (WS05). 

A useful characterisation of the temporal response of the 
medium is the intensity- weighted average delay. A, defined by 



U{t, t) = Uj 5{t — Tj) exp [2mu)jt] 



(6) 



(2) 



/„°°drT(r,t) • 

The time-average of this quantity, (A), provides us with a simple 
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Figure 5. The interstellar propagation delay, A(t) (solid line), as deter- 
mined from the holographic image shown in figure 3 with low amplitude 
coefficients (Im^I < 0.004) set to zero. Also shown is the (unweighted) 
mean delay, (A) = 15.2 (dashed line), for the observation. The dot- 
ted line shows a weighted mean (14.8 /is), with the weight for each time 
sample being equal to the intrinsic pulsar flux, /jj, as determined by the 
modelling procedure described in §2 of this paper. 

gauge of the influence of wave propagation on pulse arrival time 
for this Une-of-sight for the particular time and frequency intervals 
covered by our data. For the image U shown in figure 3 we com- 
pute (A) ~ 17 /^s. However, this value is clearly an overestimate 
because noise in the reconstruction - the noise floor is seen in the 
top panel of figure 4 - biases the estimate upward. A better estimate 
is available if we employ the image U with low-amplitude coeffi- 
cients set to zero (as per the lower panel in figure 4); this yields 
the result (A) ~ 15.2 /xs. Substantial contributions to this mean 
delay arise from the whole region < r < 120 /is, with relatively 
minor contributions from each of the discrete features - "A, B, C, 
D" of HUl et al (2005) - seen at large delay; together these features 
contribute about 10% of the measured (A). There is no contribu- 
tion from the artifacts associated with the conjugate image as these 
occur in the region t < and are excluded by the definition given 
in equation 6. 

Finally we return to the influence of the beat terms, described 
by B\ these terms cause the propagation delay A(f) to vary over the 
course of our 45 minute stretch of data. As noted earlier, the beat 
terms are not negligible and for our data we find that there are sub- 
stantial variations in A of ±40% around the mean value (A) ; these 
variations are plotted in figure 5. Large gradients in A(<;) are seen 
in figure 5 so that, for example, a pulse arrival time measurement 
at the start of our observations and one taken 5 minutes later would 
have differed by approximately 11 jis. We emphasise that the be- 
haviour seen in figure 5 is specific to this Hne-of-sight and to the 
particular set of frequencies covered in our data. We expect that the 
temporal variations in A would have been smaller if our data had 
covered a greater bandwidth than the 1.56 MHz used here. With a 
broader observing bandwidth we would have finer delay resolution 
in our image U, so we would be better able to separate components 
which have similar values of t and these separated components 
would not beat (see equation 4), so the importance of B would di- 
minish. It is beyond the scope of this paper to attempt a detailed 
analysis of delay variations; here we simply note that the type of 
behaviour seen in figure 5 is potentially problematic for precision 
pulsar timing. 



5 DISCUSSION 

Over an interval of several months we expect PSR B0834-K06 to 
exhibit changes (5(A) in the mean delay (A) as the pulsar moves 
behind different regions of the scattering screen. The electric field 
image shown in figure 3 clearly demonstrates that this particular 
screen has a very patchy distribution of scattering material so that 
observations made at other epochs are expected to exhibit quite dif- 
ferent distributions of power along the parabolic locus. As a hypo- 
thetical example we can imagine that at another epoch the holo- 
graphic image might look hke figure 3 at positive Doppler shifts, 
but show no scattered power at negative Doppler shifts; in this 
case the mean delay would be roughly half of what we have mea- 
sured. We therefore expect that over an interval of several months 
PSR B0834-I-06 will exhibit changes J(A) ~ (A). The anticipated 
delay changes (5(A) ~ 15 /xs are very large compared to the preci- 
sion with which millisecond pulsars can be timed (e.g van Straten 
et al 2001). B0834-I-06 is not a millisecond pulsar, and even if it 
were it would not be used for precision timing experiments pre- 
cisely because this Une-of-sight is known to exhibit very striking 
effects from multi-path propagation in the interstellar medium. It is, 
however, salutary to see how large the influence of the interstellar 
medium can be on pulse arrival times. Moreover the line-of-sight 
to B0834-I-06, although unusual, is not unique — B 1133-1-16 and 
B1929-I-10, for example, appear to show similar, striking effects 
(Putney and Stinebring 2005). Nor is it a very distant pulsar, so the 
scattering structures which have been revealed in the present data 
are probably abundant in the interstellar medium. B0834-I-06 is, 
however, a relatively bright source and other, fainter pulsars might 
be viewed through similar scattering media without being recog- 
nised as such because the scattered intensity is small. The data of 
Putney and Stinebring (2005) support these points as inany of the 
pulsars which they studied appear to show strong scattering arising 
from multiple, physically distinct regions, and they note that very 
sensitive observations are required in order to reveal the presence 
of these media. Furthermore the very structured nature of the scat- 
tering medium seen in figure 3 tells us that a pulsar which shows 
no measurable interstellar timing delays at one epoch could exhibit 
large delays at other epochs. Presumably interstellar scattering me- 
dia exhibit a broad range of physical properties, with a correspond- 
ingly broad range of influence on pulse arrival time measurements, 
and for any given pulsar we should not ask ourselves "whether" but 
"at what level" is an extended scattering tail present? 

A coinmonly used strategy for mitigating the influence of the 
interstellar medium on pulsar timing experiments is to undertake 
the timing observations at high radio frequencies. The rationale 
for this is based on the assumption that the scattering originates in 
distributed Kolmogorov turbulence, for which the expected prop- 
agation delay falls very rapidly with frequency. Distributed Kol- 
mogorov turbulence is not a good model for what we see in fig- 
ure 3. Currently we are not able to predict what interstellar tim- 
ing perturbation would be measured, for this pulsar and this epoch, 
at radio-frequencies outside the 1.56 MHz observing band of the 
present data. It is likely that the interstellar delays would be smaller 
at higher frequencies - because the cold plasma refractive index 
declines with frequency and so the scattering angles decrease at 
higher frequencies for any given plasma structure - but we cannot 
say how much smaller. From data spanning many months, Hem- 
berger and Stinebring (2008) have used a secondary spectrum anal- 
ysis to estimate multipath propagation delays at frequencies be- 
tween 1,150 MHz and 1,500 MHz, measuring values which range 
from one to two orders of magnitude smaller than our estimate of 
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(A). Their data refer to a different line-of-sight (PSR B1737+13) 
and cannot be directly compared with our result; however, their fre- 
quency coverage is great enough to allow them to estimate the fre- 
quency dependence of the propagation delay for their observations. 
They find that the propagation delay is consistent with a power-law 
frequency dependence of the mean propagation delay, with power- 
law index —3.6 ± 0.2 — slightly less steep than expected for dis- 
tributed Kolmogorov turbulence (power-law index —4.4). 

The possibility of large, transient interstellar propagation de- 
lays makes it prudent to quantify the interstellar delays at each 
epoch where an accurate pulsar timing measurement is desired. 
In this paper we have shown how those delays can be quantified 
when a recording of the pulsar dynamic spectrum is available. The 
technique we have presented has the merit of being able to accu- 
rately determine the relative propagation delays, even for complex 
scattering geometries, at any epoch of observation. The technique 
does have its limitations though. Foremost among these is that the 
delays are all relative measurements and the origin of coordinates 
(i.e. r = 0) must be determined by some other means; removing 
this Umitation is a high priority goal for future work. 

For cases where the great majority of the scattering arises in a 
single, thin screen we anticipate that the key to progress lies in con- 
structing the holographic image by modelling the structure of the 
phase screen, rather than forming the image U in delay-Doppler 
space as we have done in this paper. Modelling the phase screen 
is computationally much more demanding but it has several advan- 
tages over the approach we have taken here. First it allows the in- 
terstellar propagation delays to be determined on an absolute scale, 
because the propagation geometry is fixed in the model. Secondly 
it allows us to calculate the wave field for various different observer 
locations, thus permitting direct comparisons between different ob- 
serving epochs and between different telescopes (e.g. Very Long 
Baseline Interferometry) at the same epoch. Thirdly it permits di- 
rect comparison between data at different frequencies. And by the 
same token there would be no difficulty modelling data taken over 
a single, wide bandwidth even if the phase screen has large spa- 
tial variations in wave dispersion. By contrast the present approach 
is not ideal for wide bandwidth data because of smearing due to 
differential dispersive delays, across the band, from a range of dis- 
persion measures. Finally, if we model the phase screen itself the 
results immediately provide powerful constraints on models of the 
scattering medium, because the phase structure tells us the electron 
column-density structure and the Une-of-sight magnetic field if full 
polarisation information is recorded. 
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6 CONCLUSIONS 

Interstellar holography is a precise new technique which affords 
detailed insights into the interstellar propagation of radio pulsar 
signals. The holographic image reconstructed from our test data 
reveals a complex scattering structure whose physical nature is 
unknown at present. Holographic imaging can be used to deter- 
mine the influence of multi-path propagation on pulse arrival time 
measurements and thus to correct for these propagation delays. In- 
terstellar propagation delays are unpredictable and can potentially 
make large contributions to the systematic errors in pulse arrival 
time measurements. It is therefore prudent to make provision for 
holography in every instance where an accurate measurement of 
the unperturbed arrival time is desired. 
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